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The Constrained Adiabatic Trajectory Method (CATM) is reexamined as an integrator for the Schrodinger 
equation. An initial discussion places the CATM in the context of the different integrators used in the 
literature for time-independent or explicitly time-dependent Hamiltonians. The emphasis is put on adiabatic 
processes and within this adiabatic framework the interdependence between the CATM, the wave operator, 
the Floquet and the it, t') theories is presented in detail. Two points are then more particularly analysed and 
illustrated by a numerical calculation describing the i° n submitted to a laser pulse. The first point is the 
ability of the CATM to dilate the Hamiltonian spectrum and thus to make the perturbative treatment of the 
equations defining the wave function possible, possibly by using a Krylov subspace approach as a complement. 
The second point is the ability of the CATM to handle extremely complex time-dependencies, such as those 
which appear when interaction representations are used to integrate the system. 



I. INTRODUCTION 

The numerical solution of the Schrodinger equation 
ihd^/dt — H"i> is a central element in the understanding 
of experiments which involve molecular collisions or in- 
teractions between molecules and electromagnetic fields. 
For energy-resolved experiments, stationary theories are 
used. Thus, in the case of quantum diffusion theory, the 
close-coupling formalism 1 and the Lippmann-Schwinger 
approach^ solutions to the time-independent Schrodinger 
equation H^f = are sought which correspond to a 
precise total energy E and also to precise asymptotic con- 
ditions and which are found by integrating differential or 
integral equations. 

In the case of time-resolved experiments, or when 
the information on the asymptotic solutions provided 
by the close-coupling techniques is not sufficient (such 
as in laser control problems), time-dependent treatments 
are favoured. When the dynamics is driven by a time- 
independent Hamiltonian, several algorithms can be used 
to propagate the wave packet which represents the sys- 
tem in Hilbert space. These include the algorithms pre- 
sented by Leforestier et al 3 -, the second-order differencing 
scheme (SOD), the split operator method and the short 
iterative Lanczos propagation. Most of these methods 
can also be used when the Hamiltonian explicitly de- 
pends on time. However in these cases, the length of 
the integration steps must be reduced in order to han- 
dle any rapid time variations of the hamiltonian matrix. 
The propagation scheme is then based on the decompo- 
sition of the evolution operator into small increments of 
duration At: 



N-l 



U(t,Q)= Yl U((n + l)At,nAt) 



(1) 



where At = t/N and 

U(t + At,t) = exp[-(i/K)H(t + At/2)At]. (2) 

Propagation errors due to this scheme are proportional 
to (At) 3 and involve commutators of the Hamiltonian 
at successive times. These errors cancel out when the 
H matrix does not depend on time but there are also 
errors due to the approximate calculation of the action 
of exp[-(i/h)H(t + At/2) At] on the wave function ^(t). 
Thus in the three-point SOD scheme which is based on 
the equation 



^(t + At) « - At) - 2iAtH^{t)/h 



(3) 



the propagator is conditionally stable and the accumu- 
lated error per time step is equal to 



error 



(StE rn ) 3 
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(4) 
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where E m is the eigenvalue of the discretized Hamiltonian 
with the largest modulus. 

In field-matter coupling problems, the difficulty aris- 
ing from the presence of high frequencies can be circum- 
vented by introducing the Rotating Wave Approximation 
(RWA) 4 . However this approximation generates rather 
high error plateaus when the step At becomes too small, 
which is the case for intense laser fields^. Difficult com- 
bination of slow quasi- adiabatic evolutions on long time 
scales together with rapid partial or localised time varia- 
tions induce very large spectra for the hamiltonian matrix 
(e.g. in the theory of a radiative association experiment 
involving cold molecules fragments). Such cases create 
difficult problems of error accumulation for all integra- 
tors, although the amplitude and the distribution of the 
resulting errors in the spectrum is not the same for them 
all. 

To obtain high-accuracy integrators for multi- 
dimensional systems evolving adiabatically, one can also 
introduce the symplectic partitioned Runge-Kutta meth- 
ods. Using the work of Gray and Verosky on real and 
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time-independent Hamiltonians^, Sanz-Serna et Portillc 7 
have generalized the method to time-dependent Hamil- 
tonians by transforming the system into an autonomous 
equation by introducing an additional conjugate pair of 
variables (V,T = t). Another accurate time propaga- 
tion method for an explicitly time-dependent Hamilto- 
nian have been produced by Kormann 5 -, by replacing the 
Hamiltonian H in Eq.([5]) by a suitable truncation of the 
Magnus series^ H and by using the short iterative Lanc- 
zos scheme for computing the matrix-vector multiplica- 
tion exg>[-(i/K)HM\^, 

While the stability of an integrator, its accuracy and 
its ability to conserve the norm of the wave function 
are important features, it is also necessary to consider 
other elements such as the calculation time needed for 
a given accuracy, the required memory capacity, the 
complexity and generality of the integrator and also any 
constraints which could prevent its use in some cases. 
For instance some integrators such as the SOD cannot 
handle non-hermitian Hamiltonians^. The split operator 
scheme 1 " requires that the kinetic operator does not 
mix coordinates and their associated momenta. The 
multi-configuration time dependent Hartree (MCTDH) 
method^— requires important preliminary work to 
rewrite the kinetic and potential operators, while the 
calculation of higher-order terms of the Magnus devel- 
opment is a complicated task which is only tractable 
if the couplings have separated time and coordinate 
dependencies 5 . 

In this article we investigate the performances of the 
Constrained Adiabatic Trajectory Method (CATM^iS 
as a global integrator for the Schrodinger equation, with 
particular emphasis on a system which is adiabatic, in the 
sense that the system is correctly described by small sub- 
spaces spanned by eigenvectors of the molecular Hamil- 
tonian H(x,t), or by Floquet eigenstates of the field- 
dressed Hamiltonian. The CATM is well suited to the 
description of systems driven by Hamiltonians with ex- 
plicit and complicated time variations. This method does 
not have cumulative errors and the only error sources are 
the non-completeness of the finite molecular and tempo- 
ral basis sets used, and the imperfection of the time- 
dependent absorbing potential which is essential to im- 
pose the correct initial conditions. 

In sec. m the CATM theory is placed in context with 
regard to other treatments such as the Floquet theory 
or the (t,t') method, with emphasis on the concept of 
adiabaticity and on the compatibility between this con- 
cept and the time-dependent wave operator theory. Then 
three points are particularly studied, all related to the 
fact that the CATM proposes a global integrator for ex- 
plicitly time-dependent Hamiltonians and thus is not in 
the category of methods described by Eq.([2]). In sec. Hill 
we present some comparisons between the CATM, the 
SOD scheme and the split-operator method. The sec- 
ond question we ask in sec. IIVI is how the influence of a 
Krylov growing subspace algorithmic directs the conver- 
gence properties of the CATM as compared to a pertur- 



bative recursive distorted wave approximation (RDWA) 
approach to solve the wave operator equations^ (with 
or without absorbing potential, because the expansion 
of the Floquet spectrum under the influence of the ab- 
sorbing potential^ 5 - also directly acts on the convergence 
properties). In sec. El another important point emerges 
from a numerical problem that we have noted in previous 
CATM calculations^, in the case of a multistep propa- 
gation (if the time interval is too long to be treated with 
only one global step, it can be divided into several large 
steps treated in succession). Using the absorbing opera- 
tor sometimes leads to high-frequency parasites charac- 
teristic of the Gibbs phenomenon. We then try to eval- 
uate the benefits of introducing an interaction represen- 
tation before applying the CATM to implement the time 
propagation. This provides a new test for the method, 
in the presence of general time variations in the Hamil- 
tonian. Sec. IVII is devoted to the conclusion. 



II. STATIONARY AND DYNAMIC TREATMENTS FOR 
ADIABATIC PROCESSES IN MOLECULAR PHYSICS 

The main difficulty in studying adiabatic processes 
comes from the fact that in most cases adiabatic or 
even quasi-stationary molecular interactions are com- 
bined with fast time variations. A purely stationary 
process where the time only appears as a global phase 
in the wave function implies that H is self-adjoint and 
time- independent. On the other hand, as soon as the 
Hamiltonian is non-sclf-adjoint and acquires resonance 
states, characteristic times appear, such as the lifetime 
of the initial state or a characteristic tunneling passage 
timely 

Things are even more complicated when the Hamilto- 
nian becomes explicitly time-dependent, either because 
certain classical degrees of freedom are present and are 
coupled with quantum degrees of freedom or because in- 
teraction representations are used during the calculation. 
There is now no simple expression for the evolution oper- 
ator and the Dyson expansion in powers of the Hamilto- 
nian is no longer consistent with Eq. ((TJ) , except if we use 
a Magnus expansion with a time step At which depends 
on the order of the Magnus expansion. 

The difficulty of constructing the Magnus series can 
be circumvented by using the (*,*') theory^.. Thus 
the Schrodinger equation for a time-dependent Hamil- 
tonian can be solved in the same way as that for a 
time-independent Hamiltonian by working within the ex- 
tended Hilbert space K,. This extended space was first 
introduced by Sambe22 for periodic Hamiltonians and 
was then generalized by Howland^. In short, the (t, t') 
method solves the Schrodinger equation 

ihjL*(x,t) = H{x,t)*(x,t) (5) 

by adding a new variable t' to define the extended Hilbert 
space. The corresponding wave function ^(x, t',t) is re- 
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lated to V(x,t) by 

V{x,t) = y{x,t',t)\ t , =t , (6) 

with 

¥(a?, t 7 , t) = exp[-(i/»)frjr(x, 0(* - to)W(x, *o) (7) 
where i/p is a Floquet-type operator: 

H F (x,t')=H(x,t')-ih-^. (8) 

The choice between two possible representations of the 
initial state depends on the initial process studied. De- 
pending on the circumstances one can choose a time- 
independent initial state, 

V(x,t',t ) = y(x), (9) 

or an initial state which is well-defined for a specific initial 
time t , 

V(x,t',t ) = 6(t' -t )V(x). (10) 

We thus find an integration scheme which is based on 
Eqs.(JTJ) and (J2J but with a Hamiltonian which belongs 
to the larger extended Hilbert space. This can possibly 
create memory capacity problems, especially when the t' 
interval is very large. 

In the periodic case, i.e. H{t) = H(t + T) with 
T = 2tt/uj, one can transform the dynamic problem into 
an equivalent time-independent infinite-dimension eigen- 
value problem 2 ^ and generalize it to the complex quasivi- 
brational energy formalism by including finite C 2 repre- 
sentations of the molecular continual. Using the quan- 
tum variable 9 — Lot and introducing the Floquct Hamil- 
tonian: 

H F {6) = H{6)-ifk^^- (11) 
09 

which is defined in the enlarged Hilbert space K. = % ® 
C%{d6 /27r), the evolution operator acting in the enlarged 
space then becomes 

U Hp (t,t Q ) = cxp[-iH F (t - t )/h] (12) 

and is related to the evolution operator in the Hilbert 
space Ujj{t,t ) as follows: 

U HF (t,t ) = T- u tU H (t,t )T ut (13) 

where r^t = e lultd / de is a phase translation operator 
which acts on the functions of £2- We thus obtain a direct 
relation between the standard solution ^ in the Hilbert 
space H and the solution <3> in the extended space K: 

t) = Tutrix, t, 9) = t,9 + cot). (14) 

This formulation also establishes a connection between 
the quantum and the semi-classical formalisms for field- 
matter interactions at the intense field limit 2 -'. It also 



provides a way to describe the mix of adiabatic and sud- 
den effects in some experiments by using several time 
scales. 

However, the (t, t') and the Floquet theories are handi- 
capped by having quite large memory requirements in nu- 
merical applications. Moreover, in the enlarged Hilbert 
space, the calculation of the action of the operator 
exp[— iH F (t — t )/H] [Eq.(|12p] on the initial state remains 
a delicate problem when the H F spectrum is very dis- 
persed, even if the Chebyshev global scheme can be used. 

However, another approach is possible. At the adia- 
batic limit, we can stop searching for exact numerical 
solutions by adopting an adiabatic approximation such 

*(i) « exp{4 / E(t')dt'- [ (<p(t')\d<l)(t')/dt'}dt'U(t) 

(15) 

where cj>(t) is an instantaneous eigenvector: 

H(t)<p(t) = E(t)<t>(t), (16) 

and where the initial wave function ^{t ) is assumed to 
be equal to the instantaneous eigenvector <j>(t ). The 
main weak point of this approach is that it is rigor- 
ous only at the purely adiabatic limit. Such a case is 
exceptional and in most cases, the dynamics generates 
non-adiabatic couplings which mix several eigenvectors. 
Adiabatic formulae such as Eq. (fT5|) must then be gener- 
alized by introducing degenerate active spaces and non- 
abelian geometric phases. In the last part of this sec- 
tion we demonstrate that the wave-operator is undoubt- 
edly the better framework to describe this generaliza- 
tion. Indeed the non-adiabatic couplings are generated 
by the operator ihd/dt and if one renders zero this oper- 
ator in the fundamental equation which defines the time- 
dependent wave operator one obtains the fundamental 
equation which defines the stationary wave operator. The 
stationary form is then the pure adiabatic limit of the 
time-dependent form. In the following we will denote by 
==> the passage from the time-dependent to the station- 
ary equations induced by this pure adiabatic limit and 
by 7^=> the non-passage. Evidently one cannot go from 
the Schrodinger equation ([5]) to the eigenvalue equation 
(fTrJ)) by setting d/dt — > in the first one, 

(H(x, t)-ih—)V{x, t) = (H(x, t)-E(t))<f>(x, t) = 0, 

(17) 

because the operator ihd/dt also generates the dynamic 
phase which is associated with E and which is an integral 
part of the wave function. 

These drawbacks disappear within the framework of 
the wave operator theory because the time-dependent 
wave operators does not include the dynamic phase, 
which is factorized separately. We define P a as the pro- 
jector corresponding to the finite group of non-perturbed 
eigenvectors which mix under the influence of non- 
adiabatic couplings, Q being the projector onto the com- 
plementary space. Let P t be the projector associated 
with the corresponding group of perturbed eigenvectors 
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at the time t when the Hamiltonian takes the value H(t) 
and U (t, 0, H) be the evolution operator. The Bloch wave 
operator is defined b y 28 ' 29 



n f = PtiPoPtPo)- 1 =P + Q X l P 



(18) 



(The inverses are defined within the subspace S Q , 
(PoPtPo)^ 1 is the inverse of Pt within the space S ), and 
the time-dependent wave operator is defined by^ 



O(i,0) = Uit^H^PoUit^-H^o)- 1 
= P o + Q o X(t,0)P o . 



(19) 



These operators define a generalized adiabatic framework 
in which the wave function is an instantaneous linear 
combination of the eigenvectors spanning the subspace 

I 



with projector P t . The time-dependent wave operator 
factorizes the dynamic phase and the non-Abelian Berry 
phase inside U (t, 0; H e g)22-: 

U(t, 0; H)P = fl(t, 0)U(t, 0; H eS ) (20) 
with H eS (t) = P H(t)fl(t, 0)P o (21) 

These phases (which express rapid evolutions generated 
by the non-adiabatic couplings) are thus separated from 
the adiabatic evolution, the latter being included in tt. 
The result then is that the fundamental equations for 
fl(t,0) [Eq.flU])] and fl* [Eq.([l8])] adiabatically corre- 
spond to each other by cancelling out the operator d/dt. 
Instead of the non-implication (|17[) we now have: 



Q {1 - X(t, 0))H (i)(l + X(t, 0))P o = ihdX(t, 0)/dt Q a (l - X*)if(i)(l + X l )P = 



(22) 



A direct consequence is that the adiabatic limit of the 
time-dependent wave operator is given by a succession of 
instantaneous Bloch wave operators^. 

The wave operator theory is compatible with both the 
(t, t') method and the Floquet theory, and the remarkable 
property of Eq. (j2"2"|) at the adiabatic limit is conserved if 
we work in the enlarged Hilbert space JC. Thus, if H (t) is 
T = 2tt/uj periodic (T can be arbitrarily large), the two 
fundamental equations which define fl(t, 0) and f2* in the 
K space can be derived (with t as a quantum variable 
and no longer a parameter). Denoting the wave operator 
within the K, space by f2 we have the implication 



ft(H(t) - ihd/dt)Q 
= (H(t) - ihd/dt)Q 



(23) 



We note that the equations become identical by trans- 
forming Hp(t) = H(t) — ihd/dt into H(t), i.e. by ne- 
glecting the time-derivative operator (however on the left 
hand side t is a true coordinate of the K, space, while it 
is only a fixed parameter on the right hand side.) 

Nevertheless there is not a perfect equivalence between 
the two left-hand sides of Eqs.(g2]) and (gSJ). The equa- 
tion: ihdX(t,0)/dt = Q (l- X(t,0))H(t){l + X(t,0))P o 
is a non-linear evolution equation within the Hilbert 
space, which uses an imposed initial X(t = 0,0) which 
is consistent with the chosen initial wave function. Inte- 
gration of the equation indirectly gives the wave func- 
tion ^(x,t), since the initial conditions are automat- 
ically satisfied. By contrast the left part of Eq.(f23|). 
ilHp(t)H, = Hp(t)n, is a stationary eigenvalue equation 
within the extended Hilbert space and its solution gives 
Floquet eigenvectors with an intermediate normalization. 
For a degenerate S space, no column of f2 has a priori 
an initial value compatible with the initial wave function. 
For intense couplings, a large number of columns will un- 
fortunately be necessary to approach the solution of the 
Schrodinger equation (so that a subspace S with pro- 
jector P Q of high dimension will be necessary) . The left 
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equation in (1231) then provides VL in a global way, but con- 
trary to its counterpart in Eq. (|22p , which is amenable to 
a step-by-step integration, it is inapplicable to integrate 
the time-dependent Schrodinger equation. 

To overcome this major difficulty, within the frame- 
work of the CATM^ 4 a complex absorbing potential is 
added to the Hamiltonian H(t) during an extension of 
the integration interval, in such a way as to impose an 
initial vector (t = 1 compatible with the initial value 
of the wave function ^(t = 0). This can be done us- 
ing only a non-degenerate subspace S , (P )- Within a 
one dimensional subspace the wave operator represen- 
tation becomes a column. In a previous paper—, we 
have studied the matrix expressions of a general time- 
dependent absorbing potential V for any initial wave 
function >F(:z;,t = 0). In this framework, the integration 
of 

n(H F (t) + v(t))n = (H F (t) + v(t))n (24) 

gives, in a complete basis for the IC space, the expression 
of the column |f2) which is a Floquet eigenvector |A) be- 
longing to the first Brillouin zone (with an intermediate 
normalization), i.e. 



|fi) = |A)/(*,n = 0|A) 



(25) 



where \i,n = 0) is the basis vector for /C in the first 
Brillouin zone and is chosen to have maximum overlap 
with the initial wavefunction. | A) is a Floquet eigenvector 
such that 



{H F +V)|A) = E X \X) 



(26) 



where V ensures the equality between (t = 0|A) and the 
initial wave function. Eventually the wave function is 
simply given by 



|tf(t)) = exp 



-iE x t/h 



<*|A>. 



(27) 
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This formulation then gives a global solution for the 
Schrodinger equation with an explicitly time-dependent 
Hamiltonian. This generalizes finite-basis treatments 
within the Hilbert space T-L to finite-basis treatments 
within the extended Hilbert space JC. The method uses, 
in addition to the radial complex absorbing potentia l 32 ! 33 
(which is necessary to reveal resonances), extra absorb- 
ing operators asymptotically placed along the time axis 
to impose the initial conditions. Contrary to that of 
Eq. (fT5"j) where the non-adiabatic couplings are neglected, 
this solution is rigorous even if it is based on an adiabatic 
hypothesis. The addition of a time-dependent absorbing 
operator makes the wave function proportional to a single 
Floquet eigenvector [Eq. (|27]l ] which is calculated in an it- 
erative way within the framework of the wave operator 
theory. 

Similar concepts have been introduced long ago by 
Peskin et a l 34 ' 35 in the framework of the (t, t') theory, 
for the calculation of Green functions within the ex- 
tended space K. Introducing a f'-dependent absorbing 
potential to impose the boundary condition for the t' 
axis, they have shown that it was possible to replace the 
time-dependent Schrodinger equation by an inhomoge- 
neous time-independent linear system and thus to cal- 
culate transition probabilities using a scattering matrix 
formalism. They solved the linear system using a Krylov 
subspace-based iterative method in combination with a 
Fourier grid preconditioner. Despite the numerous simi- 
larities between this approach and ours, the fundamental 
working equations and the iterative procedure we use are 
different. 

From a more pragmatic point of view, the aim of the 
following sections is to give some comparisons between 
the CATM and other schemes and to analyse in detail 
how the presence of the absorbing operator affects the 
several possible variants in the method used to determine 
the eigenvector |A). 



III. COMPARISONS BETWEEN THE CATM, THE SOD 
SCHEME AND THE SPLIT-OPERATOR METHOD 

In Sections IIIII IIVI and [V] we make numerous test 
simulations with the well-known example of the pho- 
todissociation of the molecular ion H% within the Born- 
Oppenheimer approximation. This example is a simple 
ID dynamical system but it constitutes a significant test 
because the CATM is not focused, as MCTDH, on the 
treatment of multi-dimensional quantum systems with 
the use of efficient time-dependent basis-set. The CATM 
proposes a new scheme to integrate the dynamics driven 
by a complicated and fast time-evolution. These two as- 
pects (large dimensions and complicated time-evolution) 
are nevertheless sometimes correlated. The use of inter- 
mediate representations allows us to reduce the basis set 
dimensions, but in so doing it makes the time-dependence 
complicated (cf. Section [Vj). In the present section, a 
strong adiabatic laser-field envelope with several hundred 
optical oscillations is studied in order to test the capacity 



of our model to reproduce such extreme adiabatic situa- 
tions. 



A. Model for II., illuminated by intense pulses 

We study the nuclear vibration and we only take into 
account the first two effective potentials 3 -^ in the two low- 
est electronic states 2 S+ and 2 S+. The nuclear Hamil- 
tonian is the sum of two terms 



H = H + W(t). 



(28) 



Here Hq = K + Vq(x) is the field- free Hamiltonian of 
H.2 i which is pre-diagonalized on a radial grid basis us- 
ing a grid method 3 -^ in the presence of a radial complex 
absorbing potential 32 ' 33 , to obtain the vibrational eigen- 
basis with 200 eigenvalues {Ej} and bi-orthogonal eigen- 
states as well as the electric moment operator rep- 

resented by the matrix puj — We then calculate 

the dynamics in the presence of a semi-classical intense 
time-dependent electric field. In our example the electric 
field envelope is given by the gaussian function 



E(t) = E exp 



(29) 



with t = 1000 au (cf. Fig. [J). The total duration of the 
adiabatic pulse is T = 2t m = 10000 au (i.e. 0.24 ps) and 
the carrier wave angular frequency u) is 0.2958678 au, 
which corresponds to a wavelength of 154 nm. Within 
the framework of the dipole approximation, the coupling 
term W(t) is 



W{x,t) = n(x).E(t). 



(30) 



The CATM is then based on the Floquet eigenproblem 



H (x) + W(x,t) + V(x,t) - ih^- 

ot 



\X)=E X \X), (31) 



V being a time-dependent absorbing operator present 
only on the supplementary interval t 6 [To , T] and 
defined^ in Table |TJ The absorbing interval is AT ~ 
3600 au with a centred bell shape for the absorbing op- 
erator (cf. Figf2]), 



Vabs(t) = -i Vo sine 2 



(t-Kn) 

AT 



(32) 



AT 
2 



with t' m = T 

The fundamental Floquet period is the total du- 
ration T. The time description can be made us- 
ing time-periodic functions (t\n) = l/VTe~ 2 ™ f ' T 
(n e N, n = -N/2 . . . (N/2 - 1)) as a finite ba- 
sis representation (FBR), and the associated discrete 
variable representation (DVR) is defined by = 

1 l N T,n=-N/2 e- 2m(t ~ u ^ T . N is the number of Fourier 
basis functions, or equivalently the number of grid points 
which describe the time dimension. The required ex- 
tended Hilbert space can be of quite large dimension 
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Vabs(U) (column I) 

V abs (U) -^|*^.x(V r oh ,(t t )+B J -S,) 

: 

(row I) 

x (K fe (ti) + JSj - £,) Va6.(*i) 



TABLE I. Matrix representation of one block i* of the absorb- 
ing potential within the bi-orthogonal eigenbasis set {\j}} of 
H . 



if one or other component (H or C 2 ) has a large di- 
mension. The calculation of |A) can be efficiently un- 
dertaken using the wave operator theory in the case of 
a one-dimensional active space. We then have to solve 
Eq. flAll) (see Appendix), fl being simply proportional to 
the eigenvector. Eventually the transition probabilities 
— KjI^W)! 2 as wei l as the dissociation proba- 
bility 



^ iss = i- E ioi*(*))i 2 , 



(33) 



bound states 



can be calculated. 



E 
< 

CD 



03 
LU 




4000 6000 
Time (a.u.| 



10000 



FIG. 1. Adiabatic laser pulse with angular frequency uj = 
0.2958678 au and total duration 10000 au (0.24ps) with a 
gaussian shape. 
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FIG. 2. Imaginary part of the time-dependent complex ab- 
sorbing potential defined in Eq. (|32l) . 

In terms of memory, the CATM mainly requires the 
storage of a vector of the extended Hilbert space, which 



is modified iteratively. Because of the use of FFT, the 
computational effort of the program with respect to the 
time propagation scales as 2Ni ter N log(iV) where N is 
the number of Fourier basis functions and Nn er the num- 
ber of iterations needed to converge in the wave operator 
calculation (see Appendix^ . The minimum N is roughly 
proportional to the total duration for a given maximum 
frequency of the field. If N is sufficient to stabilize the 
fast Fourier transforms, then the accuracy of the results 
does not depend directly on N but rather on the absorb- 
ing potential parameters. Ni ter (10-40 in general) does 
not depend on the dimension of the matrix Hp- 



B. Comparative integration schemes 

The CATM is a global propagator for explicitly time- 
dependent Hamiltonians. However, since the Cheby- 
shev polynomial development of the evolution operator 
exp(Hp(t — to) j (ih)) in the extended space offers another 
very precise global solution^, it appeared worthwhile 
first to compare these two global solutions. In the Cheby- 
shev approach the global integrator is constructed on an 
iterative way and requires at least a minimum number 
Ncheb of iterations given by: 



Ncheb > 



AE T 

2h ' 



(34) 



where AE is the complete energy range (here the Floquet 
energy range) and T the total duration. In the present 
case, we have 



AE ~ Max(Si) + N 



2tt 



(35) 



With Max(i^) ~ 2. 934a. u. (i.e. the maximum value 
of the Ho spectrum) and N — 2048, we obtain a mini- 
mum number of iterations Ncheb — 21084. Exactly as for 
the CATM, each iteration of the Chebyshev polynomial 
scheme requires the multiplication of Hp by a vector of 
the extended space. By comparison, the CATM in combi- 
nation with the wave operator theory generally converges 
within only a few tens of iterations under the same con- 
ditions. The CATM appears then as largely competitive 
compared to the Chebyshev scheme. 

It is also interesting to compare the CATM with two 
non-global propagators, namely the SOD scheme given 
in Eq.Q and the split-operator method. To construct a 
significant comparison, the representation of the molec- 
ular H(t) for both the CATM and the SOD is made on 
the same molecular basis, the eigenbasis of Hq- Thus 
we can use the simplest form of absorbing potential for 
the CATM if the initial state is an eigenstate of Ho- Of 
course, by doing this, we lose the advantage of the more 
sparse hamiltonian matrix representation when it is rep- 
resented using a DVR grid basis for x. With the SOD, 
the accumulated error per time step is proportionnal to 
St 3 where 8t = T/N so d is the time step. Thus after N so d 
propagation steps if we wish to have a final fixed error 
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smaller than a given number e, i.e. 



N aod x I a 



T 



< e, 



(36) 



we must choose N soc i greater than (& 3 ^ 2 ^-^r 

The problem can also be studied with the split- 
operator method, which requires H (t) to be represented 
on a DVR grid basis for x. This scheme is based on a 
splitting of the kinetic and potential energies in the evo- 
lution operator for each stepi^: 



exp[ ~H8t) =exp( -~{K + V)St 



-iKSt/{2h) -iVSt/h -iKSt/(2h) 



(37) 



where V = Vq(x) + W{x,t). In the present case, a sec- 
ond splitting is necessary for the potential term because 
W(x, t), which represents the coupling terms between the 
two electronic effective potential curves, does not com- 
mute with Vq: 



exp 



i vst 



-iV St/(2h) -iWSt/h -iV 8t/(2h) 



(38) 



In the 1-D model, the kinetic energy is diagonal in 
the momentum representation (FBR), and the potential 
energy is block-diagonal in the coordinate representation 
(DVR grid for x). More precisely, exp(— iVoSt/h) is diag- 
onal in the two central blocks corresponding to the two 
surfaces. exp(— iWSt/H) possesses a diagonal representa- 
tion equal to cos(WSt/H) in the same central blocks and a 
diagonal representation equal to i sm(WSt/h) in the two 
off-diagonal blocks. Of course this property will reduce 
the split-operator computational effort with respect to 
the other two techniques, which are implemented on the 
eigenbasis of H = T + Vo(x) where the coupling W(x, t) 
is not block-diagonal (which is in a sense a more general 
exercise for the methods). 



C. Results 

We focus on the transition probabilities to the first 
bound states and on the dissociation probability. The 
system is quasi-adiabatic and mainly follows the ini- 
tial state, as is shown in Fig. [3j There are small 
non-adiabatic transitions during the pulse (for instance 



-Po->i ^ 10 ). Final inelastic probabilities are all 
smaller than 10~ 15 . 

The three methods are successively applied to this ex- 
ample with a variable number of time steps (SOD and 
split-operator) or with a variable number of Fourier ba- 
sis functions (CATM). Figj4] shows the values of the final 
dissociation probability given by the CATM and the SOD 
scheme, both working in the H eigenbasis and by the 
split-operator method working in the DVR grid basis on 
x. We have mentionned CPU-times using a small work- 
station. For an equivalent accuracy on the final Pdiss, the 
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FIG. 3. Dissociation and transition probabilities for sub- 
mitted to the pulse of Fig. [TJ The initial state was the fun- 
damental state v = 0. 



required time points numbers and CPU times are very 
different. The CATM is more efficient than the SOD 
scheme and the split-operator method is the fastest to 
converge to a correct value of the final dissociation prob- 
ability. 

Nevertheless in such a quasi-adiabatic problem it is 
also important to look at the small non-adiabatic transi- 
tion probabilities. Fig. [S] shows that the CATM and the 
SOD scheme give better results for small probabilities 
than the split-operator (which failed to correctly repro- 
duce the beginning and the end of the dynamics, even 
if a large number of steps is used). The non-significant 
decrease at the end of the process given by the CATM 
comes from an imprecision in the Floquet eigenvalue E\ 
which appears in the exponential of Eq. (I27I) . 



IV. THE INFLUENCE OF THE ABSORBING 
POTENTIAL AND OF A KRYLOV SUBSPACE 
PROCEDURE ON THE CONVERGENCE OF THE CATM 

From now on we only work with the CATM. In the 
current section we choose to treat an ultra-short pulse 
with total duration T = 212.9 a.u. (i.e. 5.15 fs), t = 40 
a.u. with the same angular frequency uj — 0.2958678 a.u. 
(cf. Fig. [6|). CATM calculations are driven with a small 
basis of 256 Fourier functions to describe the time evo- 
lution. The absorbing interval is AT = 70 au with a 
greater amplitude. 

For a given duration of the pulse, the different com- 
putational features that can be modified are studied suc- 
cessively. These features are : 

• the choice of the integration technique used to find 
the wave operator [recursive distorted wave ap- 
proximation (RDWA) perturbative calculation or 
RDWA plus Krylov subspace algorithm, see Ap- 
pendix [A] ; 
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FIG. 6. Laser pulse with angular frequency ui — 0.2958678 au 
and total duration 212 au with gaussian shape. All results in 
the current section relate to this pulse (with variable inten- 
sity). 



FIG. 4. Final dissociation probability given by the CATM 
(squares), the SOD scheme (rounds) and the split-operator 
method (triangles) vs number of Fourier basis functions 
(CATM) or number of time steps (SOD and Split-operator). 
The associated CPU-time are mentioned near each point. The 
CATM and the SOD scheme work within the Ho eigenbasis 
and the split-operator works within a DVR x-basis. 
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FIG. 5. Time evolution of the transition probability 
|(2|*(t)}| 2 given by the CATM with N = 2048 Fourier func- 
tions (cpu time: 42.5s), by the SOD method with N = 10 6 
time steps (cpu time: 471.1s), by the split-operator scheme 
with N = 10 5 time steps (cpu time: 12.9s), with N = 10 6 
time steps (cpu time: 126.7s) or with N = 10 7 (cpu time: 
1244.3s). 

• the amplitude of the time-dependent absorbing po- 
tential Vo [cf. Eq. ([32)] - If Vq is too weak the eigen- 
vector is just a Floquet eigenvector but is not con- 
nected to the correct initial condition. If Vq = 0, 
this is no longer a CATM calculation. 

• the maximum amplitude of the electric field Eq. 

Instead of approaching the Floquet eigenvector by adding 
successive correction terms to a trial vector, the Krylov 
subspace procedure uses the correction terms to con- 
struct a small but growing subspace in which direct diag- 



onalisation gives a better estimate of the Floquet eigen- 
vector. This procedure may improve the convergence 
properties of the CATM. In particular, an interesting 
point is to study how this expected improvement inter- 
acts with the expansion of the Floquet spectrum^ in 
presence of the absorbing potential; this expansion has 
an important influence on the convergence properties. 

The calculation is halted when the following conver- 
gence criterion is satisfied: 

norm [(H F - E X )\X)] < 1(T 12 . (39) 

We count the number of iterations required to reach this 
level of convergence for the two methods, which indicates 
the efficiency of the calculation and allows us to estimate 
the radius of convergence. Here we assume that the initial 
state is one of the {\j)} eigenstates, \j = i). The biggest 
(j 7^ i) component e at t — (which should ideally be 
zero) is also monitored to estimate the quality of the 
results: 

e = max [\(j\$(t = 0))| 2 ] with j ^ i. (40) 

First we determine a "non-connected" Floquet eigen- 
vector, in the absence of any time-dependent absorbing 
potential. In such a case, the calculation converges to an 
arbitrary Floquet eigenstate and not to the Schrodinger 
equation solution. Fig. [7] shows that the required num- 
ber of iterations varies roughly linearly as a function of 
Eq (with a few exceptional points) until divergence sets 
in, whatever the choosen method. The Krylov algorithm 
gives the convergence up to Eq = 0.5, whereas the simple 
RDWA diverges beyond E = 0.25. 

We then make a calculation with a quite large absorb- 
ing potential (Vq = 0.4). Fig. [5] shows a relatively linear 
behaviour for the two methods, with a clear advantage 
for the Krylov method in terms of speed but a not so 
marked one in terms of radius of convergence (Eq = 0.95 
for the simple RDWA, Eq = 1.20 for the RDWA+Krylov 
procedure). Nevertheless we must note that the qual- 
ity of the results is slightly better for the simple RDWA 
[e = 4.79 x 10 -15 when Eq = 0.3 as compared with 
e = 2.25 x 10~ 12 , see Eq.(|30"])]. Some numerical results 
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FIG. 7. Number of iterations required until convergence 
vs electric field amplitude Eq, with Vo = 0, using RDWA 
(rounds) or RDWA+Krylov procedure (triangles). All the do- 
main of convergence is covered. For each case, the last point 
(coloured in black) is the edge of the convergence domain, in 
the present calculation conditions. 
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main of convergence is covered. For each case, the last point 
(coloured in black) is the edge of the convergence domain, in 
the present calculation conditions. 



are given in Table |TT] for several runs using RDWA or 
RDWA plus Krylov algorithm for different electric field 
amplitudes to check the stability of the results. 

Fig. [9] shows more clearly the influence of the absorb- 
ing operator. We observe a divergent behaviour in the 
absence of the absorbing potential or for a weak ampli- 
tude (Vq < 0.05) when the perturbative scheme is used, 
while the Krylov scheme converges in the same condi- 
tions (to an eigenvector which is not connected to the 
initial conditions). Then, as the absorbing operator is 
introduced, the advantage for the Krylov method disap- 
pears completely and both the methods converge in a 
totally identical way. For larger amplitudes (Vo > 0.3), 
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FIG. 9. Number of iterations required until convergence 
vs absorbing potential amplitude Vo, with Eo = 0.5, using 
RDWA (rounds) or RDWA+Krylov procedure (triangles). 



the Krylov scheme again shows a significant advantage 
in terms of convergence. 
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FIG. 10. The e of Eq. (|40|l vs absorbing potential amplitude, 
with E = 0.5, using RDWA (rounds) or RDWA+Krylov pro- 
cedure (triangles). 

However, we cannot study the convergence properties 
without keeping a check on the quality of the results. As 
a measure of quality we monitor the degree of connection 
with the initial conditions^ by calculating the largest 
non-absorbed probability residue [the e of Eq. pO)) . the 
smaller e, the higher the quality]. We can see in Fig. [10] 
that the residue is almost identical for both the meth- 
ods when Vo < 0.25, but this is not a significant domain 
because the residue is too large to obtain accurately the 
solution of the Schrodinger equation. The residue regu- 
larly decreases as an exponential of the absorbing oper- 
ator amplitude, but eventually the Krylov scheme shows 
an error which stagnates at e ~ 10~ 13 while the residue 
continues to decrease for the simple RDWA. This RDWA 
advantage in terms of quality effectively cancels out the 
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TABLE II. Comparison of the CATM final transition probabilities P(\j),To) as a function of the electric field amplitude Eo 
(in units of 10 14 W.cm~ 2 ) with 2 different procedures : Simple RDWA (A), RDWA+Krylov subspace diagonalization (B). The 
biggest initial residue is defined in Eq. (|40p . The absorbing potential amplitude was Vo = 0.4. 
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previous advantage for the Krylov method in terms of 
convergence for a strong absorbing potential. 

It is also possible to compare how the two techniques 
approach the solution for a choosen calculation, with 
To = 0.25 and V = 0.4. Fig. QT] shows that the 
RDWA+Krylov procedure eventually approaches the so- 
lution faster than the simple RDWA procedure, after a 
temporary initial slowness. 

The results obtained imply that the absorbing poten- 
tial best facilitates the convergence during the search for 
the Floquet eigenvector when a simple RDWA perturba- 
tive method is used. We can explain this as follow: while 
the absorbing potential is necessary in order to constrain 
the initial condition and to obtain the solution to the 
Schrodinger equation, it also happens that under its in- 
fluence the spectrum is dilated. Under these conditions 
a perturbative approach becomes efficient, without the 
need to use variational corrections through the Krylov 
procedure. As a consequence the use of a Krylov sub- 
space algorithm for the diagonalization is not always jus- 
tified for the CATM, although in the absence of an ab- 
sorbing potential the use of a Krylov procedure gives an 
advantage in calculating a non-connected Floquet eigen- 
vector (a less important case since it does not solve the 
Schrodinger equation). When the absorbing potential is 



strong, the Krylov procedure converges faster but gives 
results of poorer quality. Since the presence of the ab- 
sorbing operator is essential in order to connect the eigen- 
vector to the initial condition in the CATM procedure, 
we conclude that the use of a Krylov subspace technique 
for the calculation of the eigenvector is not worthwhile. 



V. USE OF AN INTERACTION REPRESENTATION 
FOR THE HAMILTONIAN BEFORE APPLYING THE 
CATM 

A. The Gibbs phenomenom and the CATM 

From now on, all the CATM calculations are made 
with the RDWA procedure (cf. Appendix [A")l . Whatever 
the choice of the integration algorithm for the determina- 
tion of the eigenvector, in the general case of the absorb- 
ing potential given by Table HI some numerical difficulties 
can appear, especially in the particular case of a multi- 
step propagation (as explained in a previous article^). 
It is crucial that the terms (Ej — E{) x v f°/ v I'° which are 
present in the column no. I of the absorbing operator 
should be present only during the additional time inter- 
val [To , T] to produce the correct absorption but must be 
definitely absent during the physical time interval [0, To]. 
For this purpose we tried to use the Heaviside function 
in time, multiplying the absorbing operator by a discon- 
tinuous function equals to one on [T),T] and to zero on 
[0,T ]. In practice this produced some false or diverg- 
ing results; evidently the spectral representation and the 
numerous FFT we used are not compatible with the use 
of such discontinuous functions. Some converged results 
have been presented in a previous paper— but those re- 
sults could not be considered as completely satisfactory, 
expecially as regards the problem of the column contain- 
ing {Ej - Ei) x 

To avoid this problem, two options are possible. The 
first one is to ensure the continuity of any time-dependent 
function varying too rapidly by using smooth transition 
functions which will soften any sudden variation by pro- 
gressively turning of or turning on the sudden functions 
on the artificial interval [To > T] . This is conceivable but 
can be quite difficult to realize and may not be compati- 
ble with the aim of reproducing the exact initial condition 
exactly at the end of the interval without any interme- 
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diate continuity interval before the final time T. The 
second option which is much more simple is to work in 
the interaction representation with respect to the time- 
independent diagonal of the hamiltonian Hq (i.e. EiSu). 



B. Working in the interaction representation 

If we denote the evolution operator from t to t' asso- 
ciated with the hamiltonian X by U(t' ,t, X), then the 
evolution operator associated with H = H + W(t) can 
be written^ 

U(t, 0, H) = U(t, 0, flo) U(t, 0, V mt ) (41) 

with 

V mt (t) = U- l {t,0,H Q ) W(t) U(t,0,H o ). (42) 

Practically, we apply the CATM to the transformed 
hamitonian 



with 



Hij{t) = Wij(t) x exp I -~(Ej - Ei)t 



W ij {t) = iHj E{t), 



(43) 



(44) 



and after the propagation which gives the ^(t) we trans- 
form this result to obtain the correct wavefunction ac- 
cording to 



(45) 



01*(*)) = (j|*(*)> x ^P ( -jE,t 



Our main aim is to avoid the numerical problems previ- 
ously described, because now Ha — 0, so that we do not 
need to use the problematic column in the absorbing op- 
erator. We can also use this tranformation as a test for 
the CATM in order to see if our approach applies to sys- 
tems with various and complicated time-dependencies ev- 
erywhere in the Hamiltonian. We would also like to know 
if a treatment using an interaction representation has an 
influence on the convergence properties of the CATM. 

When we make the complete transformation of 
Eqs. (|43l) . some large terms due to the non-hermiticity 
can appear in the hamiltonian. Thus, as we work in the 
eigenbasis of Hq, some of the basis states corresponding 
to the continuum have eigenvalues with quite large neg- 
ative imaginary parts. This can be the source of numer- 
ical problems because of the factor exp (— i(Ej — Ei)t): 
if Im(Ei) ~ —1, \ e ~ l { E i~ E i) l / h \ can easily reach values 
as big as 10 80 ! However, these particular states are not 
very important in the dynamics because they are always 
completely localised at the edge of the grid where the 
radial optical potential is placed. Moreover, because of 
the Franck-Condon factors, the dipole moment coupling 
terms between such states and other states are always 



more than 1000 times smaller than the smallest cou- 
pling term between the states corresponding to eigenval- 
ues with more reasonable imaginary parts. To illustrate 
this approximation we show the spatial form of two of 
these particular eigenstates in Fig. [12] We chose to ne- 
glect these problematic transitions in this case. 
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FIG. 12. Spatial dependence of the modulae of the eigenstates 
of the field-free molecular ion : first two eigenstates of the 
first electronic potential curve 2 S^" on the left: \{x\j = 0)| 
(line) and |(z|l)| (dashed line) ; first 3 pseudo-eigenstates 
of the second electronic potential curve 2 Ej": |(a;|102)| (long 
dashes, in the middle) \{x\j = 100) | (dotted line, on the right), 
| (x 1 101) (dotted-dashed line, on the right). These two last 
functions located at the edge of the grid correspond to eigen- 
values with Im(_Ej) ~ — 1 and are some of those which play a 
very minor role in the dynamics. 

A second type of interaction representation can be en- 
visaged, with respect to only the real parts of the diago- 
nal Re(Ei), with the following separation of the hamilto- 
nian: H = Re(i/o) + ilm.(Ho) + W(t). This corresponds 
to another transformation and we apply the CATM on 
the modified hamiltonian 



Hij(t) = Wij(t) x exp I --(Re(.Ej) 
S jj (t)=iJm(E j ), 



Re(Ei))t 



(46) 



which gives the intermediate solution 'l'(i). To obtain 
the solution corresponding to the original hamiltonian, 
we must then calculate 



= (#(<)) xexp --Re(^)t 



(47) 



C. Results 

In this section the pulse is slightly longer, with a gaus- 
sian turning on and of and a plateau of 85 au. It is shown 
in Fig. [131 The other parameters are exactly the same 
as those in section llVl 

The results corresponding to computational parame- 
ters presented in Table Hill are given in Table HVl For each 
run, we can change the number of steps, we can use or 
not use the interaction representation of Eqs. (|43"| and we 
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part of H (C) [Eqs.g6 



100 150 
Time (a.u.) 

FIG. 13. Laser pulse with pulsation ui — 0.2958678 au and 
total duration 254 au (i.e. 6.14 fs) with gaussian turning on 
and off and continuous wave during 85 au. All results in the 
current section relate to this pulse with varying intensity. 

Run Inter. Repr. Number of Absorbing operator 
[Eqs.(|43|l] time steps amplitude Vb 



(1) 
(2) 
(3) 
(4) 
(5) 
(6) 



yes 
yes 
yes 
yes 
yes 



0.3 

0.3 

0.15 

0.3 

0.15 

0.3 



TABLE III. Computational parameters corresponding to the 
results of Table ITVl The electric field amplitude was 0.25 and 
the initial state is the fundamental state \i = 0). 



can modify the amplitude of the absorbing potential Vq. 
For the pulse shown in Fig. II 31 the hnal transition prob- 
abilities to some bound states |(j|\Ij(To))| 2 as well as the 
hnal dissociation probability l-£ hound ?tatcs |0'|#(T ))| 2 
are calculated for each set of computational parameters 
and the numerical results are compared. Comparison be- 
tween run (1) and run (2) confirms that the interaction 
representation is correct. Runs (3) and (4) use two steps 
of about 127 au (i.e. the time interval is divided into two 
equal parts which are successively propagated using the 
CATM) and runs (5) and (6) use four time steps of about 
76 au. Comparison between (2) and (3-6) indicates that 
the division of the calculation into several long steps us- 
ing the interaction representation is valid even if we can 
see some minor variation in the final results. In previous 
articles, the influence of the absorbing potential on the 
results has been analysed in detail. Here we can just re- 
mark that between run (3) and run (4) the increase of 
the absorbing operator amplitude increases the quality 
of the results, in the same way as it does between runs 
(5) and (6). 

Table [V] presents a comparison between the results 
given by three versions of the CATM, always working 
with only one global time step. It is associated with 
Fig. [T4J which shows the number of iterations required 
to obtain the convergence with the three different pro- 
cedures. The first curve is the simple CATM (A), the 
second uses the interaction representation with respect 
to Hq (B) [Eqs. (|4"3")l ] and the third corresponds to the 
partial interaction representation with respect to the real 
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FIG. 14. Number of iterations required until convergence vs 
electric field amplitude Eq, with Vb = 0.3, using the CATM 
and RDWA, without interaction representation (rounds) or 
using an interaction representation following Eqs. (|43|l (trian- 
gles) or Eqs. (|46[) (squares). All the domain of convergence is 
covered. For each case, the last point (coloured in black) is 
the edge of the convergence domain. 

The most simple procedure (A) shows the lowest ra- 
dius of convergence, with the first difficulties appearing 
from an amplitude of 0.35. Using the interaction rep- 
resentation extends the radius of convergence and even 
decreases the number of iteration until we obtain a good 
approximation of the eigenvector. The version (B) makes 
easier calculations in several steps involving the absorb- 
ing operator in its most general form of TableUbut needs 
more iterations than the direct calculation (A) . With the 
procedure (B) we note some small imprecisions due to 
the neglect of the large imaginary parts of the eigenval- 
ues which are located at the edge of the grid, to avoid 
exponentially increasing terms. The interaction repre- 
sentation with respect to the real parts (C) increases the 
radius of convergence without needing more iterations 
to satisfy the convergence criteria. The gain in terms 
of convergence radius and speed is remarkable with the 
procedure (C). In this one-step scheme, the gain is due 
to the fact that the procedure replaces the integral of 
the couplings by Fourier transforms of these couplings, 
at the Bohr frequencies (Ej — Ei)/h. We then obtain a 
more perturbative calculation. 

The accord between (A) and (C) is remarkable in Ta- 
ble |V] The convergence acceleration does not automati- 
cally corresponds to a gain in term of CPU time. This is 
simply because this system includes couplings n(x)E(t) 
with separate spatial and time dependencies, while this is 
no longer the case as soon as we use interaction represen- 
tations. The number of FFT required increases with the 
use of the interaction representation, but this is specifi- 
cally due to the elementary character of the studied sys- 
tem. It will not happen for more complicated systems, 
where the interaction representation will produce an ap- 
preciable gain of CPU time. 
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TABLE IV. Comparison of the final transition and dissociation probabilities to bound states with the different conditions 
described in Table UTTl 



Electric Field Procedure Biggest initial Final Pdiss 
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TABLE V. Comparison of the CATM results as a 
function of the electric field amplitude (1 corresponds 
to 10 14 W.cm- 2 ) with 3 different procedures : Simple 
CATM (A), CATM+interaction representation [Eqs. (|43l ] 
(B), CATM+interaction representation with respect to real 
parts [Eqs. (|46[ )] (C). The biggest initial residue e is defined in 
Ea. ([40"jl . The absorbing potential amplitude is Vo = 0.3. 



VI. CONCLUSION 

The present article proposes a global integrator for ex- 
plicitly time-dependent Hamiltonians. The addition of a 
time-dependent absorbing potential in order to impose 
the initial conditions in the Floquet treatment has the 
added advantage of improving the convergence in the so- 
lution of the Floquet eigenvalue equations, mainly be- 
cause of the expansion of the spectrum. 

The introduction of a Krylov procedure, which turns 
out to be very efficient for the calculation of Floquet 
eigenvectors without complex absorbing potential, only 
offers a minor advantage when the absorbing potential 
is present. Moreover a Krylov approach needs to store 
many vectors within the extended Hilbert space, which 
can be a problem with the CATM because the spaces 
in question have quite large dimensions. The absorbing 
potential then provides an advantage equivalent to the 
benefit of a Krylov subspace procedure but with a stor- 
age which is limited to one vector. 

The second point of interest is that the CATM is con- 
sistent with the use of interaction representations. In 
this framework, the global CATM integration reduces the 
non-diagonal coupling amplitudes (for a relatively slow 
perturbation). In general, for step by step integrators, 
the presence of supplementary oscillations due to the in- 



teraction representation would be a supplementary diffi- 
culty, but with the global CATM scheme this facilitates 
the convergence. 

We insist on the interest of the global approach 
to integrate systems which show complicated time- 
dependencies, possibly due to the use of an interaction 
representation. For half an oscillation, 4 grid points are 
sufficient in FFT calculation to obtain accurate results 40 , 
whereas step by step integrators can require hundreds of 
points, this number becoming larger if the required prob- 
abilities are small. Nevertheless, until now our test sim- 
ulations on H% or other small systems with the CATM 
often take approximately the same CPU time compared 
to optimized step by step integrators, and depending on 
the conditions the CATM is sometimes faster. But our 
algorithm is still in work and the CPU time optimization 
is only one of the numerous motivation to continue the 
study of the CATM. Indeed this algorithm keeps its ad- 
vantages already described in previous articles 1 ^"—, espe- 
cially the ability to the repetition of calculations, or the 
complete control of the accuracy of the results depend- 
ing of the characteristics of the absorbing potential and 
of the Fourier basis set. Moreover the global structure of 
the CATM seems to be compatible with parallel compu- 
tation, which could significantly increase the interest of 
this algorithm. 
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Appendix A: Determination of the time-dependent wave 
operator 

To obtain the results presented in this article, the 
constrained Floquet state has been calculated using the 
wave operator method. A basic RDWA iteration proce- 
dure was first applied, and then we tested the combina- 
tion of this procedure with a Krylov subspace construc- 
tion with approximate diagonalization. These techniques 
have some similarities with Davidson's method 41 -. 
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RDWA iterative method 



with 



The fundamental equation for the wave operator f2 for 
the Hamiltonian Hp is^ 



fl-W = (i _ X (n) )H F {l + X (n >) 



(Al) 



For the one-dimensional case, fi is directly proportional 
to one eigenvector of Hp. For the general case, has the 
form Q = P + X with X = Q XP , Po and Qo being 
respectively the projectors on the active space So and on 
the complementary space Sq. Finding f2 is then equiva- 
lent to finding X. The effective Hamiltonian is defined by 
H c s = PqHfU. From Eq. (|Al|) . it is possible to get a self- 
consistent equation for X, with H' an arbitrary diagonal 
matrix in the complementary space (H' = QqH'Qq) and 
|/) a vector of the finite basis representation such that 
Qo|/> = ]/>: 

(f\XP Q = ((/| [H F -H']X + (f\H F P ) (A2) 

x (P H eS P - (flH'lfiPo)- 1 . 

One possible choice (namely the Recursive Distorted 
Wave Approximation) is 

H' = (Q (l - X)H F (1 + X)Q Q ) diag (A3) 

= (Q (l - X)H F Q ) diag . 

This choice leads to the equation 

(f\XP = (A4) 
(f\H F XP - (f\(l - X)H F \f)(f\XP Q + (f\H F P 
(H eS -{f\(l-X)H F \f))P 

which can be solved using several numerical procedures. 

2. Generalized Newton- Raphson Procedure 

We thus have to solve an algebraic equation X = 
F(X), X being a vector and F the operator on the right- 
hand side of Eq. (|A4[) . We define a trial vector l' ' (for 
example equal to the initial wavefunction delocalised over 
all the time interval) and then modify it iteratively by the 
addition of small quantities, 



X( n ) = x^" -1 ' + AX^ 



(A5) 



until we obtain a sufficiently accurate satisfaction of 
Eq.flTnj, when calculated with n( n ) = P + X^ n \ If 
the trial vector and the diagonal matrix H' are well 
chosen- 2,43 , following a classical linear procedure, the 
small increment is calculated as 



AX (n) = F(X 



(n-l)s 



(A6) 



This procedure converges almost linearly. In the partic- 
ular case of a one dimensional space Sq, Pq = |cv)(a|, we 
obtain an approximate Newton-Raphson procedure: 



(a\HV)\a) - (f\HW\f)' 
(A7) 



(A8) 



The result is directly tested and if the convergence crite- 
rion is satisfied we stop the calculation. 



3. Krylov-type procedure 

Another point of view can be adopted^: the proce- 
dure of Eq. (|A5p defines a subspace of growing dimension 
which is spanned by the sequence of vectors X^ and 
which contains progressively more information about the 
solution X. The growing Krylov-type basis set {|ej)} is 
constructed by Gram-Schmidt orthonormalization of the 
sequence of the correction terms Al'"': 



|e )=X(°) 
|ui) = AX« 



k-l 



XeolAXW, | ei) 



\/{ui\u\) 



(A9) 



)=AX^-J2h)(e i \AX^, \e k ) = 



\u k ) 



\/{Uk\Uk) 



After a few iterations, we assume that the Krylov sub- 
space almost contains the required eigenvector and so 
wish to recombine the generated approximations into 
something better. After k iterations, the Krylov subspace 
K, k C C fe is of dimension k and is spanned by the orthogo- 
nal basis Vfe G C nxfc (containing the |e») in columns). We 
can then calculate a good approximation of the solution 
by diagonalizing the restriction of Hp to the subspace 



i.e. 



H Ftk Y = EY 

with Hp t k 
and Y, E (diagonal) € 



(A10) 



V T k H F V k € 



f^kxk 
kxk 



The T superscript and the bar denotes the transpose op- 
eration and the complex conjugation. On returning to 
the original Hilbcrt space, we then obtain the following 
approximate expression for k possible eigenvectors asso- 
ciated to the eigenvalues E (N denotes the dimension of 
the extended Hilbert space): 



z k = V k Y e 



iNxk 



(All) 



Among these k vectors, we just have to identify the 
"good" one X<- k \ which gives the maximum value of the 
overlap integral of its projection at t = with the initial 
given wavefunction. This result is in general expected to 
be closer to the exact solution X than the direct iterative 
result. 
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